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A hard-sphere (HS) Bose gas in a trap is investigated at finite temperatures in the weakly- 
interacting regime and its thermodynamic properties are evaluated using the static fluc- 
tuation approximation (SFA). The energies are calculated with a second-quantized many- 
body Hamiltonian and a harmonic oscillator wave function. The specific heat capacity, 
internal energy, pressure, entropy and the Bose-Einstein (BE) occupation number of the 
system are determined as functions of temperature and for various values of interaction 
strength and number of particles. It is found that the number of particles plays a more 
profound role in the determination of the thermodynamic properties of the system than 
the HS diameter characterizing the interaction, that the critical temperature drops with 
the increase of the repulsion between the bosons, and that the fiuctuations in the energy 
are much smaller than the energy itself in the weakly-interacting regime. 
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1. Introduction 

The thermodynamic properties of a trapped interacting Bose gas have been in- 
vestigated extensively in the last decade HI El El H] [5l [6]^ g^^j. gjj^QQ Bosc-Einstcin 
condensation (BEG) was realized for the first time in a trap in 1995 Finite-size 
effects ^ and the possibility of BEG in one and two dimensions ^ have been ex- 
plored. Numerical 1^ studies of interacting BEGs in traps as well as quasi one- and 
two-dimensional condensate expansion due to the interactions have been carried 
out. Inparticular, the effect of the interaction on the transition temperature Tc 
[TT] [T^l H] has been probed extensively. Further, the noninteracting trapped Bose gas 
has been a fascinating research topic HH El H] [16] g^j-^^ thermodynamic prop- 
erties such as the BE occupation function, specific heat capacity, internal energy, 
and entropy have been calculated. A basic issue in this context is the role of the en- 
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ergy fluctuations. For example, GiorginiH^ found that, below these fluctuations 
cause a shift in the excitation frequencies of order of a few percent. 

The aim of this paper is to undertake a full investigation of the thermodynamic 
properties for a trapped hard-sphere (HS) Bose gas with small and large numbers 
of particles by using a discrete-states approach. The role of the energy fluctuations 
in determining these properties will be underlined. 

We use the static fluctuation approximation (SFA)llSllSI2Q].^^jjj(,j^ includes in a 
relatively simple manner the above energy fluctuations. The principal feature of SFA 
is the replacement of the square of the local energy fluctuation operator (AE)^ with 
its mean value ((A£^)^). [((Ai?)^)]^/^ turns out to be much less than the chemical 
potential ^ which, according to Giorgini places the fluctuations in o ur system 
in the collective-mode category. SFA has been used to study liquid "'He^^I]^ li quid 
^He^^, spin-polarized atomic hydrogen ^I^, spin-polarized "^He-He II mixtures ^i^, 
ferroelectrics and nuclear matter All these are homogeneous systems. In 
contrast, SFA is used here for the flrst time to investigate an inhomogeneous-density 
system: a trapped Bose gas whose inhomogeneiety is a result of external trapping. 
The non-uniform Bose condensed gas has been investigated by Ohberg and Stenholm 
^ in the framework of a Hartrce-Fock (HF) formulation. However, HF theory fails 
at temperatures near the critical temperature Tc, whereas SFA works well below 
Tc- Hence our present SFA formulation of the trapped Bose gas may be regarded as 
a complementary method to that of Ohberg and Stenholm. Further, an approach 
is invoked according to which the particles are distributed in discrete energy levels; 
whereas most previous work, e.g. ^ ^ ^^Sl^ used a continuous density of states. 

We find that thermal fluctuations, for the ground state of a BEG in a harmonic 
trap, reach a maximum at the BEG transition temperature, and decline afterwards 
as the temperature T rises. In addition, we find that a.t g < 1 x lO^'^ the energy 
fluctuations are very small ^--^0(10"^) and do not influence the thermodynamic 
properties very much. The entropy is found to increase with T up to Tc, beyond 
which it begins to stabilize, reaching a plateau at higher T. Hence, the highest order 
in the system occurs when T — > 0; whereas at higher T disorder is independent of T. 
It turns out that the critical temperature for the onset of BEG decreases with the 
increase of the repulsion between the bosons. It is further found that a change in the 
number of particles of the system has quite a strong influence on the thermodynamic 
properties. The Bose-Einstein (BE) occupation function is studied as a function of 
T and the effect of interaction on the distribution of the bosons in different HO 
states is investigated. We find that indeed a large number of discrete HO states is 
needed 100) to describe strongly-repulsive BEGs. 

The paper is organized as follows. In Secl2]the method is presented briefly. In Sec. 
[3] the performance of an updated version of the SFA code is analyzed. In Sec. |4]our 
results are discussed and compared to previous literature. In Sec. [5] our conclusions 
are listed. Appendix I A. 1 1 explains the SFA iterative procedure used, and Appendix 
lA. 21 describes the SFA method. 
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2. Method 

We consider an interacting hard-sphere (HS) Bose gas of N particles in an 
isotropic, spherically symmetric, harmonic trap of trapping frequency w^o = 27rx 10 
Hz and at finite temperature T. The strength of the interatomic interaction is char- 
acterized by the s-wave scattering length as, which is equivalent to the HS boson 
diameter in the low-energy and long- wavelength approximation. We consider a hard 
core repulsive interaction between the bosons. In what follows, we present only the 
key points of the method; details are relegated to Appendices lA.ll and IA.2I 



2.1. SFA 

Here we shed light on the small parameter of our approximation. As stated in the 
Introduction, we replace the square of the local energy fluctuation operator {AE„i)'^ 
for each state m with its mean value {AEm)^ ■ The energy fluctuation AEm is given 

by 

AE,n = E,n — {E,n), (1) 

where Em is the energy operator of state m [Ea. (jA.ll[) ]. and {Em) is the aver- 
age value [Eq.(|10|)]. We define the small parameter of the approximation, written 
ipF(jn,T), which appears in the thermodynamic functions in Sec l2. 51 below, as 

^%{m,T) = {{AEmf). (2) 

That is, 

{(AEmf) = {El) {Em)\ (3) 

as usual. Essentially, i^i^(m,r) is proportional to the magnitude of the interaction 
constant g as well as the average correlations between pairs of number fiuctuations 
[see Ea.(|A.4p below]. 



2.2. Hamiltonian 

To evaluate the energy of the system, we use the following mean-field (MF) 
approach. We begin with the general form of the Hamiltonian: 



H = 



2m ^ ^ ' 



V'(r) + - / drip'' {r)iP^r)Vintir ~ r')i]j{r)^l,{r). 



(4) 



To describe the interaction between the bosons in the trap, we first consider a 
two-body HS contact potential given by 



Vi„t{ri - r2) = gSiji - r2). 



(5) 
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where g = Airh^as/u is the interaction parameter, with the s-wave scattering 
length in the low-energy and long-wavelength approximation and u the reduced 
bosonic mass. We further write the field operators as linear combinations of creation 
and annihilation operators: 

■^W = ^4>m{^)hm, (6) 
m 

(j)n (r) being single-particle noninteracting wavefunctions for a particle at position r 
and bm is the boson annihilation operator. 

Substituting ([5]) and dH) into (|4]) and taking for 0„(r) the harmonic oscillator 
wavefunctions in Cartesian coordinates: 

(/)„i(r) = 0„Ja;)0,„Jy)(?!)„^(z), (7) 

with 

^ra{x) = \ ^ CXp(-xV2)i/™(x), (8) 

2™m!-y/7r 

where Hmix) is the Hermite function of order m (and similarly for y and z), one 
gets for the second-quantized Hamiltonian for interacting bosons in a harmonic 
trap 



H = ^fta;,,o(m + 3/2)6]„S„ + 

f E ^l,&l.^c3&c.y'c(r)C(r)'/'c3(r)0c.(r)dr, (9) 

C1C2C3C4 

m and q being integers representing harmonic oscillator (HO) states. Then, as 
shown in Appendix lA. 21 [from Eq. (jA.7|) to (|A.12[) ]. the MF average energy per state 
m can be written: 

^ 00 

{E^) = huJho{ni + i/2) + -g^c{n,m){hk), (10) 

n=0 

where {hk) is the BE occupation function: 

("fc> = '^^ \m — 7' (11) 

cyiY>[{Ek - /^)/3J - 1 

with /3 = l/ksT, kg being Boltzmann's constant, a statistical weight defined 
by the degeneracy of the system, and c(fc, m) the elements of the interaction matrix 
obtained from the integral in Eq.®: 



c{n,m) = -J- 



E 



n 



exp(-2u^)ij2^ ("O^m. (wj)du» 



(12) 
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Here m must fuUfil the condition m ~ nix + + "m-z , which leads to a degeneracy 
in each state, and dm = {fn + l){m + 2)/2 is the degeneracy of state m: nii = m,x, 
771,2 = TOy, ^TT-s = and similarly for n; ui ~ x, U2 = y, M3 = z- The reason c(n, m) 
has been divided by the degeneracy dm is that Eq. (|12l) is a sum over all possible sets 
of degenerate HO states {rux, my,mz) satisfying m — rrix + my + and interacting 
with all possible sets of degenerate states (nx,ny,nz) with n = Ux + Uy + Uz- 
However, Eq. (fTO|) describes a particle in a given state m for a single set {mx,my, niz) 
interacting with all the other sets (n^, ny,nz). Therefore, to ensure correct counting, 
c{n,m) must be divided by dm- Bn-.mi stands for the normalization factor for each 
direction i, and is given by 

Bni.nii — r>„ .4_,v, . i i ■ (1*^) 

Our computational limitations have restricted the rank of the matrix c(n,m), i.e., 
the maximum values of n and m. 

We use trap units for the energy and length, hiuho and aho ~ ^/fi/muho, respec- 
tively. Therefore, 

9 -> TT— = 47ra,aL; (14) 

and similarly E — > E/huho', by multiplying 5 by a factor aj^^ arising from the triple 
integral of the interaction matrix, g becomes Anas in trap units (ag/aho — >■ ««)■ 



2.3. Numerics 

The integral s ill Eq. (|12p are determined by applying a standard Gaussian quadra- 
ture tech nique The Hermitc functions are evaluated by using the recursion for- 
mula 1221 

Hn+i{x) = 2xHn{x) - 2niJ„_i(a;); (15) 

with Hq{x) = 1 and Hi{x) — 2a;, one can iterate the above formula to get Hn{x) 
for any n. 



2.4. Chemical potential 

In what follows a new, simple technique is presented that facilitates a speedy and 
accurate evaluation of the chemical potential fi for each temperature T. Since the 
sum of {fik) [Eq. pip ] for a given /i yields the total number of particles, we define 
the function 

I exp[((£^m) - n)f3] - 1 J 

which has been minimized at each T so as to tuneu(r) to a desired number of 
particles N by using Powell's minimization routine Here M is the total number 
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of states used. However, for the evaluation of fi in SFA, the 'corrected' occupation 
function riQ{'m,T) weighted by dm, which contains the fluctuations ipp [Eq. (|A.2ip ]. 
has been used instead of (n™) [Eq. (fTTj) ]. The evaluation of (pp using the SFA iter- 
ative procedure is explained in Appendix lA.ll 

Thus, one first evaluates fj, for the noninteracting system by setting {E„i) = 
htOhoi'm + 3/2) in Eq. p6|) : next this fj, is used in the evaluation of (Em) for an 
interacting system via Eq. (fTO|) . The energy of this interacting system is then fed 
back into Eq. ()16p to get fi for the intera ctin g system. This technique is somewhat 
similar to that used by Napolitano et al. M is tuned further to get the correct 
value for the specific heat capacity in the classical limit (Sec. 14.6^ . 



2.5. Thermodynamic functions 



In this section the equations used in the SFA program are reviewed briefly. The 
pressure is evaluated in trap units using 



hWho 

where 17 is the volume of the system, Q is the Gibbs partition function, and 

M 

InQ = - ^ qo{m,T), 



(17) 



(18) 



qoim, T) being a discrete function given by Ea. (jA.37p . 

According to Rochinl^, the "volume" Q. in the case of a non-uniform trapped 
Bose gas system can be written Q, = '-^ho', where D is the dimensionality of the 
system (in our case D = 3), and louo the trapping frequency. Further, we note that 
instead of evaluating P only, we have actually evaluated the 'pressure' as P x O 
which has units of energy. Thus, even if P does not have the units of force per unit 
area, just as SI does not have the units of volume, the most important fact is that 
P X Q. has units of energy. In addition, we show in Appendix IA.3I that our Eq. p8)) 
yields the thermodynamic potential for a non-uniform trapped Bose gas system in 
the limit of an infinite number of states and very small fluctuations. 

The internal energy is given by 



U 



dlnQ 
dp 



fJ,,ipF,E 



M 

E 

7n—0 



dqo{m,T) 
dp 



(19) 



where the derivative (|T9|) has been calculated at constant fi, (fp, and E. Otherwise, 
if we included the derivatives of /i, ipp, and E with respect to /?, we find that a 
physical behavior for the thermodynamic properties is not obtainable. 
The specific heat is evaluated using 



df 



N 



M 

dT ^ 

m=0 



dp ■ 



(20) 
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Finally, the entropy is given by 

S ^ "-±ip^^ ,21, 

In evaluating S, we find that fiN constributes only a small amount to S and can 
therefore be ignored. 

The thermodynamic equations (|17p and (1211) in our paper were previously used 
for nonuniform gases by Sevincli and Tanatar bu t in different units than ours. 
Further, Grether et aL^used Ea.([T7| and Rochin^^HlTj^sgfj (j2l|) for nonuniform Bose 
gases. 



3. Analysis of the performance of the SFA code 

In what follows, a brief analysis of the perfor mance of a more advanced version 
of the SFA program than that used inEB E2] [23] jg presented. This version has been 
written in the C programming language and new features have been included that 
have enabled us to obtain /Lt(T) as in Sec l2.41 as well as ipp(m,T), and perform a 
convergence check of the iteration procedure. 



3.1. Determining the Chemical Potential 

Figure [U shows /i as a function of T/Ti"' for various g, where Ti°' = 4.5113 
nK has been chosen as our reference temperature. These values of n{T) have been 
obtained by the technique of Sec l2. 41 using ■r]o{m,T). Tc°^ has been evaluated from 
the noninteracting result 

Tj°) ^ 0.94-^iVi/3, (22) 
Kb 

for N = 1000 and ujho = 27r x 10 Hz; this value is used throughout the rest of the 
paper. Below the critical temperature, /i is — as expected — equal to the ground- 
state energy of the harmonic oscillator in the weakly-interacting regime, namely, 
^ 1.5 in trap units. Above Tc°\ /i decreases steadily and acquires negative values. 
The chemical potential seems insensitive to the changes in the values of g indicated 
in the figure. 

For comparison purposes, IJ,{T) is also displayed for g = 1 x 10"'', as obtained by 
the MF technique in Sec. l2.4l using (nfc(T)). SFA and MF results shown are identical. 
The fact that /i below ri°'' 1.5) is close to the ground-state energy indicates that 
most of the particles reside in the condensate at T < Ti"-* . These results are in line 
with those of Ohberg and Stcnholm ^ and Ketterle and van Druten ^ for the ideal 
Bose gas. However, when g is increased to order 10^^, below ri"'' attains values 
higher than 1.5. 



3.2. T-dependence of ^pp 

Figured] displays the fluctuations ipp(rn ~ 0,T) for the ground state (m ~ 0) of 
systems with g = 1 x 10^^, 2 x 10^^, and 5 x 10^^. The fluctuations increase with 
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Fig. 1. The chemical potential /i for a trapped HS Bose gas of A'^ = 1000 particles as a function 
of T for different g using SFA and MF techniques. The number of states is M = 45, and the 
trapping frequency is Lijf^^ = 27r X 10 Hz. The dashed line represents fi ~ 1.5 and Ti°^ is the critical 
temperature for a noninteracting Bose gas. 
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Fig. 2. SFA fluctuations in the energy ipp{m = 0, T) as functions of T for a system of N = 1000 
trapped HS bosons and M — ofof.io o+ fk^.^.^ jr^f .^^.^^f « '-pi°) 
a noninteracting Bose gas. 



: 45 states at three interactions g. ' is the critical temperature for 



T and show a maximum at the transition temperature Tc, beyond which they are 
damped. Note that refers to the transition temperature in the interacting case; 
one can infer from Fig. [2] that Tc < Tc°^ . 
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Fig. 3. SFA fluctuations Lpp{m = 0, T) as functions of the number of iterations Niter for the 
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Fig. 4. SFA fluctuations ipp(m = 0,T) as functions of the number of iterations Niter for the 
ground state of a trapped HS Bose gas oi N = 1000 particles and M = 45 states at g = 1 x 10"* 
and different T. 



3.3. Convergence of SFA 

To check the convergence of the iteration procedure outhned in Appendix lA.ll 
</Ji?(m,T) is plotted versus the number of iterations Nuer in Figs. [31 |31 and [S] for 
several T and g. Figures [3] and IH with g = Ix 10~^ and 1 x 10""', respectively, show 
the need for a large number of iterations (~ 25) for SFA to converge to a steady 
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Fig. 5. SFA fluctuations <fip{m = 0,T) as functions of the number of iterations Nuer for the 
ground state of a trapped HS Bose gas of N = 1000 particles and M = 45 states at g = 1 X 10~^ 
and different T. 
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Critical temperatures as functions of N for a HS Bosc gas in a trap with g = 1 X 10 



(Tc ; open circles) and for a noninteracting HS Bose gas (T^ 
and MF approaches, respectively. 



solid circles), evaluated within SFA 



value of LpF{m = 0,T). In Fig.[5l with g = 1 x 10"^, on the other hand, Lppi'm.T) 
converges rapidly after only 4 iterations for all the temperatures indicated. Thus, 
the higher the interaction, the larger the number of iterations needed to achieve 
convergence. It is, therefore, of the utmost importance to check the convergence of 
SFA calculations. 
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3.4. Critical temperature 

Figure [6] shows a comparison between the critical temperature for an interacting 
(using SFA) and a noninteracting (using MF) system as a function of N. The sohd 
circles pertain to Tc°^ of Eq. (j22p for a noninteracting, and the open circles to Tc 
for an interacting, system with g = 1 x 10^"'. The SFA results have been obtained 
from the peak positions of the specific heat capacities. 

A considerable difference is noted between the two results, which shows that the 
critical temperature drops with increasing g. Further, the trend in Tc is almost the 
same as Tc°\ The critical temperature Tc obtained by SFA increases with N, as 
expected on physical grounds. 

3.5. Failure of SFA at higher g 

It is foimd that SFA breaks down at g > 1 x 10^'^. The reason for this is that, at 
the higher interactions, the fluctuations become very large. As a result the exponent 
{Em) ~ l-i- ^ becomes negative in 



This causes the second term on the right-hand side of rjo{m, T) to acquire a negative 
sign as well. If a large number of states is negative, the sum of 770 (to, T), weighted 
by the degeneracies dm, may also turn out to be negative; whereas it is supposed 
to yield the total number of particles for an optimal Essentially this kills the 
calculation, since f{pL,T) [Eq. (|16|) ] will be far away from the true minimum. We 
may call this the 'sign problem' in the SFA method. 

4. Results and Discussion 

In this section the interrelationship between N, M, and g, the energy as a function 
of T, the occupancy of states, and the thermodynamic properties are presented. 

4.1. Interrelationship between N, M, and g 

Figure [7] displays the dependence of N on M for two values of g. For each N, the 
value of M is adjusted to the classical limit of the specific heat capacity, Cv — SNks, 
in the weakly-interacting regime. One can see that, for a given g, an increase in N 
necessitates a rise in M; although the rise in M is fairly small for each increment of 
N, namely 100. Generally speaking, then, one would not need an infinite number 
of states to accommodate numbers of particles of order 100. It is when one seeks 
systems with N ^ 10^ that the number of states is expected to shoot up to very 
large values. It is further found that the change of M with N is independent of g 
in the weakly-interacting regime, and M is mainly determined by N. 



770 (to, T) 



2 1 exp[^((i^„i) - fi +'fim)] - 1 exp[/3((^m) - /i - (^,„)] - 1 




(23) 
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Fig. 7. Interrelationship between the number of particles N and the number of states M of 
trapped HS Bosc gases at two values of g. The points (TV, M) pertain to the classical limit 3Nkg 
of the specific heat capacity. 
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Fig. 8. SFA ground state fluctuation ipp{m = 0, T) vs g at several T for = 1000 particles. 

4.2. Energy fluctuations 

Figure [8] displays (pp{m ~ 0, T) versus g for iV = 1000 at various T. LpF{m = 0, T) 
rises almost linearly as g is increased. Further, if pirn ~ 0, T) rises with T up to the 
transition temperature beyond which it drops. 



4.3. Energy 

Figure [9] displays the MF energy per particle {Em=o)/N in the condensate state, 
as obtained from Eq. (fT0|l . for = 1000 particles and various g. The energy rises 



I 
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Fig. 9. The MF energy per particle for m = as a function of T for A'^ = 1000 trapped HS 
bosons, M = 45 states and three interaction parameters g. 



initially at a faster rate below Tc than above it. As expected, the energy also rises 
with g. However, even though g and the changes in it are small, the variations 
in vs T are relatively considerable; for example, at T/ri"^ = 1.5, one observes 
AE > 0.05. This indicates that a small change in g causes {nk{T)) to attain a 
noticeably stronger T-dependence, i.e., {iik{T)) is very sensitive to variations in g. 
As the temperature rises, more particles are thermally excited out of the condensate 
to higher HO states. As a result, the energy arising from the interaction of the 
condensate with higher excited states increases with T , causing an increase in the 
total energy. 

Figure [TU] displays the behavior of {Ern=o)/N when g is fixed and N is varied. 
Clearly, the energy rises with A^; as N increases, more particles occupy higher HO 
states. Notice that a relatively considerable change in N is required in order to 
increase E vs. T. For example, one needs AA'' = 300 to get AE = 5 x 10~^ at 
r/ri°' = 1.5. This shows that a chan ge in N has a weaker effect on the energy than 
a change in g. 



4.4. Comparison to analytical results 

The total mean-field energy of this work is 

AI 

Etotal = (24) 

m=0 

where (Em) is given by Eq. ()10p . Figure [11] displays a comparison between the MF 
Etotai/N (Eq.Ell) for N = 100, 5 = 1 x lO'^, and M = 100, and the analytical 
Etotal /N as obtained by Eqs.(47) and (67) of Su et al.^^ioi the same parameters N 
and g. The two results almost match for T < Tc°\ and good agreement is obtained 



October 19, 2010 0:15 WSPC/INSTRUCTION FILE sfa"paper'14Jun09 
14 





Fig. 11. The total MF energy per particle E-tatal/^ for an interacting trapped HS Bose gas of 
N = 100 particles, M = 100 states, and g = 1 X 10"'^. The open triangles show the analytical 
results by Su et al. and the solid circles the MF results of the present work. 



for T > Tr. Figure [12] is the same as Fig. [Til but for N = 1000. Again, the 
two results almost match for T < Tc°^; but beyond that, the deviation from the 
analytical results is substantial, particularly above T = 1.5Tc°\ The reason for 
this large deviation goes back to the limited maximum number of states that we 
can use. Had we been able to apply M > 100, the MF result in Fig[T2| would have 
approached the analytical result as M was increased since the excited particles 
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would have been able to occupy higher states. In Fig.[TTl the number of states 
was sufficient to accommodate all the excitations, and this is why there is good 
agreement above Tc°\ Further, the good agreement below Tc°^ arises from the fact 
that, at T < Tc°\ all particles are in the (m = 0) condensate state. Thus, there is 
no need for a large number of states to describe the system below Tc°\ In other 
words, the effect of excitations becomes significant only beyond T = Tc°\ 

Thus, the MF model developed in Secl2] works well. Likewise, the minimization 
technique of Sec l2.4l for evaluating the chemical potential is quite effective. A super- 
computer should enable one to reach states much higher than M = 100, particularly 
for the 3D case. 



4.5. Occupancy of states 

Figure fT3l shows the fractional occupation number (condensate fraction) {hQ)/N 
versus T of bosons in the HO ground state, evaluated using SFA [crosses, 770 — 
0, T), Eq. (|AT2T|) ] and MF [open circles, (nfc=o), Eq. ^] for a system of iV = 1000 
particles, M — 50 states, and g = Ix 10~*. The SFA riQ{m = 0, T) is identical to the 
MF {nm=o{T)) because the fluctuations in the energy tpp arc much smaller than the 
energy itself ai g ~ 1 x 10^"*. One observes that, slightly above Tc°\ {ho)/N drops 
to zero. In Fig. [14] we display for the same system the fractional occupation number 
for the first excited HO state {fiij/N versus T. One can see that, while {no)/N 
decreases from 1 to almost zero at Tc^\ {hi)/N increases steadily to a maximum 
at T 0.9Tc^^ above which it drops systematically to zero. This kind of behavior 
has also been observed by Kctterlc et al. Going further, we present in Fig. [15] 
("■50) /N for the same system. The occupation number in this rather high state keeps 
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Fig. 13. SFA (crosses) and MF (open circles) fractional occupation numbers {•ho/N) as functions 
of T for the HO ground state, m = 0, of a system of N = 1000 trapped HS bosons, A4 = 50 states, 
with g = 1 X 10"''. 



rising with T. If M — 50 had not been the highest state assigned in this particular 
calculation, (n^o) would have reached a maximum beyond which it would have 
declined to zero, thus allowing other particles to be excited to states higher than 
m = 50 at higher T. Hence, in general, the maximum of is expected to shift to 
higher T as m becomes higher. It turns out, then, that M chiefly controls the overall 
distribution of N particles among the discrete states. It is to be emphasized that 
when one restricts M (to, e.g., 50), the BE occupation number at any given state m 
along with the chemical potential will be overestimated. This is because the bosons 
are shuffled back to lower states and are redistributed among them, particularly for 
g > 0.01. In this case, any states higher than m = 50 remain practically unoccupied. 

Figure [T4l shows that, when T reaches Ti°\ about 2.2% of the particles occupy 
the first excited state; and FiglTS] shows that only ^ 0.5% occupy the state m = 50. 
But as T rises further above Tc^\ the particles leave the condensate state m = 
and occupy very high states such as to = 50; for example, at T /Tc^^ = 1.8 we see 
that (nso) /N ^ 0.02, i.e., 2% of the particles occupy our highest state. Nevertheless, 
for the values of g used in this investigation (5 < 1 x 10"'^), M = 50 seems to be a 
reasonable value, since {h^o)/N is only 2% in Fig.lTSl 

Further, the maximum values of N and g we are able to use, namely N = 1000 
and g = 1 X 10"'^, correspond to Nas = Ng/{A-K) = 7.96 x 10~^, which is much 
less than 1. Consequently, we are way below the Thomas- Fermi regime defined by 
Nas > 1. 
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Fig. 14. SFA (crosses) and MF (open circles) fractional occupation numbers {fii/N) as functions 
of T for the HO ground state, m = 1, of a system of N = 1000 trapped HS bosons, A4 = 50 states, 
with 3 = 1 X 10~*. 
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Fig. 15. SFA (crosses) and MF (open circles) fractional occupation numbers (h^o/N) as functions 
of T for the HO ground state, m = 50, of a system of N = 1000 trapped HS bosons, M = 50 
states, with 3 = 1 X 10~*. 



4.6. Thermodynamic properties 

4.6.1. Critical Temperature 

Figure [TBI displays {no)/N for a fixed N = 1000 particles and different g. Tc does 
not change appreciably between g = 1 x 10~^ and 1 x 10~^; but on increasing g 
to 0.05, Tc drops noticeably by —0.2Tc°^ which is about 10 times the result of 
Giorgini et al.^. AT 0.016Ti°\ the shift in Ti°^ being given by these authors 
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Fig. 16. The MF fractional occupation number {ho)/N for A'^ = 1000 trapped HS bosons and 
M = 100 states, for g = 0.00 (noninteracting system) and g = 0.05 (interacting system). 



as 



= _1.3^iVi/'^. (25) 



Thus, there is a large discrepancy between the shift of Tc as obtained from our 
MF approach and that from Eq. ((25|) . The reason for this large discrepancy is that, 
at higher interactions, the current HO wave function —originally designed for a 
weakly-interacting system— becomes inadequate for the description of a trapped 
Bose gas with g = 0.05. This is because the bosonic wave functions begin to 
broaden substantially away from their width described by exp(— in Eq.®, 
and one should rather apply a parameterized Gaussian with a variable width de- 
scribed by exp(— ax^), where a is an adjustable parameter obtainable by variational 
techniques. We thus anticipate that, if an optimized wav e function is used in our 
calculations, Tc will approach that of Giorgini et al. We will investigate this 
issue in the future. 



4.6.2. Internal energy 

Figure [17] displays the internal energy per particle U /N as & function of T for dif- 
ferent g, keeping = 1000 fixed. U rises with T; but it does not change appreciably 
with the variations in g. This indicates that the rate of two-body collisions is low 
in the weakly-interacting regime. The rise of U with T is partly due to the thermal 
expansion of the Bose gas towards the edges of the external trapping potential. As 
it approaches the edges of the trap, the gas gains more potential energy. Figure ITSl 
shows the variations of U{T) with N , g being kept fixed at 1 x 10~*. A change in 
TV causes a much more pronounced change in the thermodynamic properties than 
a change in g within the displayed range. This is because, as rises, the density 
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Fig. 17. The SFA internal energy per particle U /N as a function of T for TV = 1000 trapped HS 
bosons, M = 45 states, and different g. 




of the system rises, and therefore the rate of cohisions increases, ahering the ther- 
modynamic properties substantiaUy. In this case, U{T) is initiahy insensitive to the 
changes in g and N at lower T up to T - 0.75ri°^; but then towards T ^ T^°^ the 
deviations in the behavior of U{T) become noticeable at various TV. 



4.6.3. Specific heat 

Figure fT9l shows the specific heat capacity per particle Cy/N versus T in units 
of ks for a system of = 1000 bosons and different g. The behavior of Cy = 
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Fig. 19. The SFA specific heat capacity per particle C^/N in units of kg for = 1000 trapped 
HS bosons for different g and M. The shift in the transition temperature is indicated. 



{dU/dT)Q conforms with the behavior of U for the same systems in Fig.[T7]at g = 
1 X 10~^ and 1 x 10~^. Initiahy, the slope of U increases with T up to Tc, after which 
it begins to dechne and then stabihzes at higher T towards the classical regime. 
The appearance of a peak in can also be directly attributed to the behavior 
of fJ,{T) in the neighborhood of Tc displayed in FigH] since given by Eq. ipO)) is 
connected to /x(T) via qo{m,T) [Eq. ((X37| ]. Note that below T^, d^/dT - 0, and 
when it approaches Ti"-* , d^/dT suddenly drops to a negative value (see Fig. [1]) . This 
abrupt change in /i(r) and dfi/dT causes a jump in qo^m^T) and dqo{m^T)/dT 
which explains the appearance of a peak in Cv Alternatively, this 'bump' in Cy 
is clearly an 'order-disorder' transition since = T(dS/dT)n. We note that C„ 
stabilizes at approximately the classical value of '--^ 3NkB for g = 1 x 10^^ 
and 1 X 10~^. However, for 5 = 1 x 10"'^, a larger number of states M = 100 
was taken to accommodate the need of particles to occupy higher HO states. This 
leads to a larger classical limit for C^, namely, QNks, as we are now approaching 
the strongly-interacting regime. In fact, it was found that, for large interactions 
g >1 X 10"'^, SFA calculations do not converge for M < 100. Figure [T9l shows also 
that the signature of the transition —i.e., the peak of — shifts to lower T by an 
amount AT/t]"-* ^ —0.0186 as g rises from 1 x 10"'' to 1 x 10"^. Using formula 
(|25l) . we get for AT/T^°^ between g = AttUc = 1 x 10"^ and 1 x 10"^ the value 

Ar/ri°^ 2.944 X lO"'' K which is two orders of magnitude smaller than our 

SFA resuh for AT. 

Figure [20l shows for a fixed g and various N. We observe that the amplitude 
of the peak in C„ rises with increasing N. Further, the peak shifts to higher T as 
N is increased, indicating that Tc incr ease s. This is in line with the result for dilute 
Bose gases t]"^ - OM{hLUho/kB)N^^^ 
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Fig. 20. The SFA specific heat capacity per particle Cv/N in units of kg at a fixed g = 1 X 10 * 
and different TV and M of trapped HS Bosc gases. 



4.6.4. Entropy 

Figure [21] displays the entropy per particle S/N in units of fcs versus T for 
different 5 at a fixed N = 1000. The entropy rises steeply below Tc^\ but then 
kinks over towards a plateau. Figure [22] is the same as Fig. [ST] except that we keep 
g fixed at 1 x 10~^ and vary N. Thus, as for C^, S eventually reaches the classical 
limit at higher T. Again the change in g from 1 x 10~'^ to 1 x 10~^ does not have 
an appreciable effect on the thermodynamic properties. The entropy reveals higher 
order as T — > 0; in fact, as T — >■ ri°^ the previously mentioned order-disorder 
transition is observed. 



4.6.5. Pressure 

Figure [23] displays the pressure P{xn) as a function of T for the same systems 
of Fig. [21] P rises with T steeply up to Tc^\ after which it kinks slightly towards a 
lower slope. Figure[24]is the same as Fig. [23] but for fixed g and various N. We note 
that, at lower T up to T/T^°^ ~ 0.75, P is the same for all N; but as one approaches 
Tc^'' the pressures for different N begin to deviate substantially from each other. At 
higher T, the slope of Pn vs. T in Fig.[23l evaluated between T/ri°^ = 1.0 and 1.5, 
is ^ 2.305 x lO^'^ K~^. This value is very close to the classical value NkB divided by 
the trap energy unit hojho, where NkB/fujJho = 2.084 x 10^^ for ujho = 27r x 10 
Hz and = 1000 particles. The same is true for Fig. [24] where the slope is different 
for each A^ in the classical limit. 
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Fig. 21. The SFA entropy per particle S/N in units of fcs as a function of T for A'' = 1000 
trapped HS bosons, M = 45 states, and different g. 
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Fig. 22. The SFA entropy as a function of T in units of fc^ for a fixed g = 1 X 10 * and different 
N and M of trapped HS Bose gases. 



4.7. Comparison to previous work 

Ohberg and Stcnholm ^ presented qualitatively a numerical study of an inhomo- 
geneous Bose gas by solving the Hartree-Fock equations for systems with repulsive 
interactions and restricting the number of particles, much as we do in our present 
study. They calculated the condensate density profiles and the thermally excited 
parts of the density as functions of temperature, particle number, and interaction 
strength and found that, above the critical temperature, there is no condensate. 
They further noted that HF theory docs not describe the Bose gas well for all T 
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Fig. 23. The SFA pressure (P X fl) as a function of T for Af = 1000 trapped HS bosons, M = 45 
states, and different g. 
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Fig. 24. The SFA pressure (P X Q) as a function of T for a fixed g = 1 X 10 * and different A'' 
and M of trapped HS Bose gases. 



and that their theory breaks down near the critical temperature, contrary to SFA. 
In fact, we have shown that SFA functions well below Tc- The HF equations were 
solved by using an iterative approach which involves a density n(r), then solving 
them for this density and using the new eigenfunctions to calculate a new den- 
sity, and so on. This procedure was then iterated until one found a self-consistent 
solution for the density and eigenfunctions. The convergence of this method was 
checked by tracing the value of the ground-state energy. Thus, their main achieve- 
ment was to obtain the density of the g function of T. We note that, as we 
do, Ohbcrg and Stcnholm used a limited number of states which they did not show 
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explicitly. In contrast to their approach, wc used a fixed eigenfunction and obtained 
the thermodynamic properties by iterating the energy fluctuations [Eq. (jA.33[) ]. 

Ketterle and van Druten ^ investigated BECs in traps for a finite number of 
particles in an ideal Bosc gas. They found that the signature of BECs. namely a 
discontinuity in the specific heat capacity for a small number of particles, is very 
similar to an infinite number of particles; that as one tightens the confinement, 
Tc rises; and that BEC is possible in ID and 2D. Further, the population of the 
first excited state for oo is negligible, even in the absence of interactions. 

As mentioned is Sec. 14.51 our occupation function {ni)/N agrees with that shown 
by these authors in Fig. 2 of their paper; but they did not display the occupation 
function for states higher than m — 1. In contrast, we have displayed the occupation 
function for m = 50 in Fig. 1151 

Giorgini et al. presented results for the temperature-dependence of the con- 
densate fraction and the dependence of the critical temperature on the interaction 
between the bosons in a harmonic potential. They found that the thermal depletion 
of the condensate is enhanced in the presence of repulsive interactions, and that 
Tc shifts to lower values with increasing repulsive interactions. In addition, finite- 
size effects do affect the thermodynamic properties; however, interactions play a 
more profound role as they lead to significant depletion of the condensate. Further, 
when one increases the number of atoms, the interactions become more profound. 
Interestingly, comparing our figures displaying the thermodynamic properties for 
different A^ and fixed g to those for different g and fixed A^, we arrive at the same 
conclusion. Finally, Giorgini et al. argued that the shift in the critical temperature is 
the combined effect of a finite number of particles plus the interactions. They found 
that the shift 5Tc°^ arising from finite-size effects is always negative and vanishes 
at large A^. 

Kao and Jiang 2^ provided a theoretical treatment for the effect of interactions 
on Tc- They noted that a finite number of particles usually drags a system away 
from the thermodynamic limit and that the effect of interactions is the hardest 
to treat among all other effects. This is substantiated by the present work since, 
as mentioned in Sec. 13.51 SFA breaks down at 17 > 1 x 10~^. In particular, they 
referred to the work of Giorgini et al. El and noted that these authors incorporated 
only thermal effects into their calculations of the shift of Tc and ignored the effect 
of the condensate component. Accordingly, Kao and Jiang considered this effect 
and underlined the importance of the interaction between the condensate atoms 
themselves, since even near the transition a small fraction of atoms resides in the 
condensate state at the trap center. Inspecting our Fig. [131 we also see that at Tc 
there is indeed a remaining small condensate fraction. The interaction between this 
remaining fraction and the thermal component is considered explicitly in our work. 
Further, they treated the system of condensate plus thermal component as a two- 
fiuid model and showed that the energy gap between the thermal and condensate 
components dete rmines the shift in the transition temperature. 

Kao et al. ^ studied ideal Bose gases trapped in a generic power-law potential. 
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using the grand canonical ensemble. Corrections to the thermodynamic properties 
due to a finite number of particles were evaluated, such as the transition temperature 
of the system as well as the condensate fraction. It was further argued that when 
the system has a finite number of particles, it is very important to consider the 
occupation fraction of excited states. 

Napolitano et al. considered a finite number of noninteracting bosons trapped 
in a three-dimensional isotropic harmonic oscillator trap, and evaluated the corre- 
sponding heat capacity. For this purpose, they used the approach of Ketterle and 
van Druten In order to provide a signal for a phase transition — a discontinuous 
heat capacity — they calculated the heat capacity numerically for a fixed number of 
particles. Cat, which at the time their paper was published had not been measured. 
First, they evaluated /^(T) from the normalization condition 

oo 

TV = ^7„,?(i?„), (26) 

where r]{En) is the BE occupation function, 7„ = (n-|-l)(n-|-2)/2 is the degeneracy 
of the HO state n, and En = nhajho is the energy eigenvalue. For T < T°, the 
ground state 7i = of the system is macroscopically occupied. Once the chemical 
potential had been evaluated, the other thermodynamic properties followed easily. 
To obtain /i(T), they found the roots of the equation 

Q 

Si^i,T) = Y,lnViEn) ~ N, (27) 

for which S'(/i, T) = 0; a technique which bears some resemblance to ours in Sec. 
12.41 In our case, we adjust the value of M (their Q) to give the correct classical 
limit of the specific heat capacity for noninteracting trapped Bose gases, namely, 
SNks- However, in their approach, they ensured the convergence of fi(T) by it- 
erating Ea. ((27)) . with S{^,T) set to 0, each time increasing the number of states 
Q, calculating the resulting fJ.{T) thereof, and comparing it to that for the previ- 
ous values of Q. If the new fi{T) differed from the previous value by only a small 
increment, then this procedure was likely to have converged. Thus, they kept on 
increasing Q until convergence was achieved. Then they evaluated the specific heat 
capacity for several values of N and showed that there was a discontinuity in the 
thermodynamic limit; whereas the discontinuity vanished in systems with a finite 
number of particles. Our results for the specific heat capacity for finite systems in 
the weakly-interacting regime also reveal the absence of a discontinuity in Cy. In 
fact, the work of Napolitano et al. is one of the rare investigations of trapped Bose 
gases that use a discrete-states approach. 

Grossmann and Holthaus SI looked into the BEG of relatively few particles in traps 
and considered particularly the effects of particle number on the condensation tem- 
perature and specific heat capacity. They evaluated the thermodynamic properties 
and found that the emergence of the A— point in their specific heat capacity was 
solely due to the external trapping potential. 
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Gnanapragasam et al. ^El calculated the occupation number of atoms in the 
ground state of an ideal Bose gas and an interacting Bose gas in a trap. From 
this, they investigated the effects of temperature and interaction on the conden- 
sate properties. They found that the depletion of the condensate rises with the rise 
in interaction and temperature, and that repulsive interactions stabilize the above 
systems; whereas in the case of attractive interactions, the systems up to certain 
negative g are stabilized by the zero-point repulsive energy arising from the con- 
finement. We may argue that the signature for the stability of our systems is the 
convergence of SFA calculations to stable values of the fluctuations (pF{m,T) and 
the chemical potential n{T). In any case, we did not have to worry about the sta- 
bility issues of our systems. This is because the repulsive HS interactions between 
the bosons balance the compressional forces arising from the external trapping po- 
tential that try to squeeze the Bose gas radially towards the center of the trap. 
Further, we found similarly that as the temperature and interaction are increased, 
more and more particles are depleted out of the condensate. This result is in line 
with Gnanapragasam et al. ^i^. 

Grether et al. ^ described the HO trapping of bosons and fermions in lower di- 
mensions, the goal being to understand these structures. Further, they considered a 
d-dimensional noninteracting boson or fcrmion gas trapped by several mutually per- 
pendicular harmonic oscillator potentials and showed that, in the thermodynamic 
limit, the specific heat capacity for a trapped Bose gas in 3D reaches 3NkB. They 
found that BEG can only occur when d + S > 2, where d is the dimensionality of the 
system, and 6 the number of HO potentials used that are mutually perpendicular to 
each other. In summary, the authors calculated the thermodynamic properties and 
densities for ideal Bose and Fermi gases trapped in 5 mutually perpendicular trap- 
ping potentials. In turns out that the behavior of some thermodynamic properties 
obtained here by using SFA is quite similar to theirs (Fig. 2 of their paper). 

5. Conclusions 

In summary, we have presented a formal investigation of the thermodynamic 
properties of a trapped, interacting HS Bose gas of finite size in the weakly- 
interacting regime, using the static fluctuation approximation (SFA). Our basic 
goal was to investigate the effect of energy fluctuations on these thermodynamic 
properties. 

We find that these fluctuations are much smaller than the energy itself in the 
weakly-interacting (g < 10~^), low-density regime. When the interactions are of 
order g = 10~^, the energy fluctuations rise substantially and SFA breaks down in 
this highly-interacting regime. 

Further, the critical temperature Ti*^^ drops with increasing interactions for 
a fixed number of particles and rises with increasing N at fixed interaction 
strength. 

In addition, changing g in the range from to 1 x 10^"^ does not have any 
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significant effect on the thermodynamic properties; but changing N while keeping 
g fixed has a more profound effect. 

Wc have also found that the specific heat capacity for a fixed g increases with 
increasing N. The specific heat in the weakly-interacting regime, with a correctly 
adjusted total number of states, reaches the classical limit SNks at sufficiently high 
temperatures. 

Appendix A. 

A.l. Iterative procedure 

Here we outline briefly the iterative method used to determine the energy fluc- 
tuations ipF(jn,T) using SFA equations ^I^, modified for a trapped HS Bose gas 
system. The basic point is to solve the nonlinear equations presented below, par- 
ticularly Eq. (|A.4p . so as to obtain stable values for (pp{m,T). The procedure is 
constructed from a loop over all temperatures T, with chosen steps AT, and an 
inner loop, conducting iterations at each T to determine ippijn^T). Hence, for each 
T: 

(1) one first assigns initial values for Lpp(rn,T) for all states m and the true 
correlations between the number fluctuations (An,„Anfe)c for m ^ k. Then 
an inner loop is started which conducts an iterative procedure to optimize 
the values oi Lppim^T) and (Afi„iAnfc)c. 

(2) At the beginning of each iteration, the function f{iJ,,T) [Eq. (fT6|) ] is mini- 
mized with respect to the chemical potential /x(T) only. 

(3) Then having, determined fJ.(T) from the minimization, one evaluates the 
mean of the square of the number fluctuations: 



where 771 (to, T) is given by (jA.2ip below, and c(m, fc) is the interaction 
matrix, Eq. (fT2|) . 

(4) The correlations between the fluctuations are updated, using ((An„j)^) in 



(5) (AhjnAhk) is used to update the true correlations between the fluctuations. 




(A.l 



(A.2) 
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according to 

M 

{AnjnAnk)c = rji{m,T) g^c{m,i) {AfiiAfik) 

1=1 

= rii{m,T) gc{m,k) {{Ahkf) + 
M 

T]i{m,T)g^c{m,i){AniAnk)c- (A.3) 

(6) The fluctuations of the energy arc calculated using the updated correlations 
{AhmAnk)c- 

M 

[</5_F(m,T)]2 = —^-—y^c{m,i){Ah^An„^)c (A.4) 
771 (m,r) f-^ 

(7) One goes back to step (2), using the updated Lpp{m,T) and {AhjnAnk)^, 
and steps (2)-(7) are repeated for a specified number of iterations. 

(8) At the end of the iterations, the procedure changes the temperature and is 
repeated. 



A.2. SFA 



In this appendix, we derive the thermodynamic properties using the static fiuctu- 
ation approximation. We begin with the Heisenberg representation of the creation 
operator &]„(r) given by 



blir) = exp{TH)bl{0) exp(~riJ). 
This satisfies the equation of motion 



dr 



(A.5) 



(A.6) 



If we decompose the Hamiltonian H in ^ into a noninteracting and interacting 
part, Hq and Hint, respectively, where Hq is the first, and Hint the second term on 
the right- hand- side of Eq. ^ , then Hq commutes with b^^ according to 



H„,bl ^ najho{m + 3/2)bl 



for the interacting part, we get 



Hin t , ^i-n 



1 



-g ^ c{ni,n2,n3,m)bl^bl^b 



"3 I 



(A.7) 



(A.8 



where 



3(711,712, 713, m) = J dr(t>l^{r) (t)l^{r) (pnsir) (t>m{r). (A.9) 
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The energy Em can then be evaluated from the commutation relation 



H, hi 



— Em.- 



which leads to 

E,n = hhJho{m + i/2) + ^ c(ni, m, ns, m) 6jj^6„3 . 

"1113 

The average energy {Em) becomes, then, 

{Em) - 

fiLJhoim + 3/2) + -g ^ C(ni,m,n3,m)(6j,^&„3)(5„i„3, 



which leads to Ea. ([TU|) . 

In SFA, we define the energy operator by 

Em = {Em) + ^Em, 



(A.IO) 
(A.ll) 



(A.12) 



(A.13) 



where AEm is the fluctuation operator of the energy. The basic idea of SFA is to 
replace the square of the fluctuation operator with its mean value: 

(AEm)' « ii^Em)') = vl- (A.14) 

Going back to Eq. (jA.6p . one can write the solution to this equation as 

^L(^) = &L(0) exp[(^„.)T] cxp(A4„t). (A.15) 

With the aid of the identity 

B[a + bAEm] = r]n{m) + r]i{m) AEm, (A.16) 

where 

Voim) = ^{B[a + bipm] + B[a - bfm]} ; 



'71 (™) = {B[a + bLpm] - B[a-bipm]}; 

2(Pm 



(A.17) 



we can rewrite Eq. (jA.15|) as follows: 
blir) = 

bliO) cM{Em)T) 



cosh((p,„r) + s'mh{tpmT) 



(A.18) 
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Using the identity 

= {BA) = iTr[cxp(-/3i?)Bi], (A.19) 

where Q is the Gibbs partition fu ncti on, and replacing /3 with time r, we can obtain 
the so-called long-range equation 

(ri^i) ^ 770 (m) (i) + 771 (m) (A^,„i), (A.20) 

where, by invoking the chemical potential via H ^ H — 



Vo[m) = - < -—7— 7 + 



?7i(m) 



2 [cxp[l3{{E.m) - fl +ip,n)] - 1 exp[l3{{E,n) - H - (p,n)] - 1 
1 f 1 1 



2 (/3m [ exp [/?((£;,„) - ^ +</?,„)] - 1 CXp[l3{{E,n) - l-l - fm)] - ij 

(A.21) 



The long-range equation is now applied to derive expressions for the fluctuations 
in the number of particles and the correlations between them. Substituting A = 1 
into this equation, and using the fact that the quadratic fluctuations are symmetric, 
i.e., {AE„i) ~ 0, we obtain for the particle distribution 

(n™) = 770 (m). (A.22) 

We also define the fluctuations in the number of particles An„i via 

Aflrn = flm - {flrn), (A.23) 

which then gives 

(An™i) = 771(771) (A^™i). (A.24) 

Substituting A ~ An^ into the long-range equation, we arrive at the expression 

{AhmAhk)c = 771(771) (Ai^„Anfc), (A.25) 

where the index c denotes true correlations, i.e., m ^ k. 

Now the fluctuations in the energy arise from the interactions between the par- 
ticles. It is straightforward to show that 

AE^ = ci^p, m)Anp. (A.26) 

p 

The true correlations between the fluctuations can be decomposed into 

{AhpAhk) = ((A%)2) Sp^k + {AhpAhk)c- (A.27) 
To find ((A77fe)2), we invoke Eq. (|A.19p with A = bl to get 

{bUl3)hnk) = {hfikbl). (A.28) 
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Applying the commutation relation 



hh. we find that 



mmnk) = (1 + K + 2hk). 



(A.29) 



Since the energy operator Em commutes with h\, and Sfc, Eq. (jA.28P can be written 
in the form 



(A.30) 



Setting = fik on the right side of (|A.30p . and equating it to the right side of 
(|A.29p . we get 



1 + 2nk 



cxp[l3{Ek - m)] - 1 

(nfc)(l + 2{hk)) - 2 77i(fc) (hkAEk). 



(A.31) 



From Eg. ljA.Bip one can determine the fluctuations in the number of particles: 



{{An,)') = {nl) {nu? 

{{Afik)') = {fik) (1 + {fik)) + 2?7i(fc) L ^c(fc,p) {AfipAfik), 



(A.32) 



Setting A ^ AE in Eq. (|X24)) . we obtain an expression for the fluctuations in the 
energy, cp„: 



771 (to) (^^ = {An„iAE„i) = y^c(fc, m)(An„,A 



(A.33) 



Next, the Gibbs partition function is derived which contains the fluctuations. 
This is given by 



Q = Tr{exph/3(i7-7VM)]} 



E 



exp 



-/3Y.iEp-fi)i 



II E '^^P ~ 



p rip 



(A.34) 



where N = ^ iVp Taking the logarithm for convenience, we get 
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In Q = In < ^ exp -/3 {Ep - ^)np 

V 

-El- 



1 — cxp 



-f^{Ep - m) 



(A.35) 



Substituting E^^i = + ^^^-rm a-nd using the identity (f06|) . we arrive at 



InQ = [(jo(p,r) + qi{p,T)^Ep 



(A.36) 



Considering the symmetry of the two eigenvahies of the energy fluctuation operator, 
we finaUy get 



Ing = -^qob,T) 



1 — exp 



(3{{Ep) -11 + ^p)] } 
P{{Ep) - ^p)] }) . 



(A.37) 



From qq (p, T) all other thermodynamic properties follow. 

A. 3. InQ in the limit M — >■ oo 

Assuming, then, that (j)p{m^T) — > 0, we write first 



ln( 



M 

LQ«-^ln{l - exp[-/3(^,„(r)-/i(r))]}. (A.38) 

m=0 

Following Grether et al. 1^, we evaluate the above sum by letting M — > oo and using 
their Eq.(6): 



[exp(/3(/.(r) - ^hw/2))f 



(A.39) 



where 



= fesTlnjl - exp[-/3(3ntj/2-^)]} 



(A.40) 
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is the ground-state contribution to the grand potential Here we have used BE 
statistics and harmonic trapping in 3D. Using their Eq.(7), we obtain 



with a = (/i — 3huj/2)/ (kBT). Apart from the fact that the treatment of Grcther 
et aZ. Pluses a reference potential f^o and a reference chemical potential 3huj/2, this 
thermodynamic potential is the same as Eq.(l) in the paper of Rochinl2Hl-\\r];ic]-c the 
'harmonic' pressure is also given by P = —fl/V, with V = ui^^ —the 'harmonic' 
volume. 

In conclusion, then, Eq.(5) in our treatment with (f)F{m,T) — >■ and M —>■ oo 
gives the thermodynamic potential for a non-uniform trapped Bose gas. Hence all 
the thermal properties that follow from it are also suited for non-uniform systems. 
But as stated before, we were only able to use M ~ 100 because of our computational 
limitations. 
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